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Abstract 

We solve the homogeneous Bethe-Salpeter (HBS) equation for the scalar, 
pseudoscalar, vector, and axial-vector bound states of quark and anti-quark in 
large Nf QCD with the improved ladder approximation in the Landau gauge. 
The quark mass function in the HBS equation is obtained from the Schwinger- 
Dyson (SD) equation in the same approximation for consistency with the chiral 
symmetry. Amazingly, due to the fact that the two-loop running coupling of 
large Nf QCD is explicitly written in terms of an analytic function, large Nf 
QCD turns out to be the first example in which the SD equation can be solved 
in the complex plane and hence the HBS equation directly in the time-like 
region. We find that approaching the chiral phase transition point from the 
broken phase, the scalar, vector, and axial-vector meson masses vanish to zero 
with the same scaling behavior, all degenerate with the massless pseudoscalar 
meson. This may suggest a new type of manifestation of the chiral symmetry 
restoration in large Nf QCD. 
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1 Introduction 



Spontaneous chiral symmetry breaking is one of the most important properties to un- 
derstand the low-energy phenomena of QCD in the real world. This chiral symmetry 
is expected to be restored in QCD at several extreme conditions such as QCD with 
a large number of massless quarks, large Nf QCD (see, e.g., Refs. [111211311100111) 
and QCD in hot and/or dense matter (see, e.g., Ref. jl]). In Ref. j2], based on 
the infrared (IR) fixed point existing at a two-loop beta function for a large num- 
ber of massless quarks (Nf ^ y^c) |TJ, it was found through the improved ladder 
Schwinger-Dyson (SD) equation that chiral symmetry restoration takes place for Nf 
such that Nf* < N f < ^iV c , where Nf' ~ AN C (= 12 for N c = 3). Then, in Ref. 
this chiral restoration at Nf 1 was further identified with "conformal phase transition" 
which was characterized by the essential singularity scaling. Moreover, such chiral 
restoration is also observed by other various methods such as lattice simulation 0, 
dispersion relation 6J, instanton calculus jjj, effective field theoretical approach 0, 
etc.. 

More attention has been paid to the property of the phase transition. Especially, 
it is interesting to ask what are the light degrees of freedom near the phase transition 
point in the large Nf QCD: For example, in the manifestation of the chiral symmetry 
restoration a la the linear sigma model, the scalar bound state becomes a chiral 
partner of the pseudoscalar bound state and becomes massless at the phase transition 
point. On the contrary, in the vector manifestation (VM) [TU1 ITT] obtained by the 
effective field theoretical approach based on the hidden local symmetry model |12J, it is 
the vector bound state which becomes massless as a chiral partner of the pseudoscalar 
bound state. Besides, from the viewpoint of the conformal phase transition j^j, it is 
natural to suppose that all the existing bound states become massless near the phase 
transition point when approached from the broken phase (see Ref. [T5]). 

Then, it is quite interesting to study which types of the bound states actually 
exist near the phase transition point, and investigate the critical behavior of their 
masses directly from QCD. Such studies from the first principle will give us a clue to 
understand the nature of the chiral phase transition in large Nf QCD. 

A powerful tool to study the bound states of quark and anti-quark directly from 
QCD is the homogeneous Bethe-Salpeter (HBS) equation in the (improved) ladder 
approximation (see, for reviews, Refs. fiJEIlEi])- When the mass of the quark is 
regarded as a constant, we can easily solve the HBS equation by using a so-called 
fictitious eigenvalue method [T^j. However, for consistency with the chiral symmetry, 
the quark propagator in the HBS equation must be obtained by solving the SD 
equation with the same kernel as that used in the HBS equation jTTl UHl l2~Uj . 
and as a result, the quark mass becomes a certain momentum dependent function. 
Then, in order to obtain the masses and the wave functions of the bound states, it is 
necessary to solve the HBS equation and the SD equation simultaneously. 

When we try to solve these two equations in real-life (Nf = 3) QCD, however, 
we encounter difficulties. First of all, for the consistency of the solution of the SD 
equation with QCD in a high energy region, we need to use the running coupling which 
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obeys the evolution determined from QCD /3-function in the high energy region (see, 
for reviews, Refs. HE]). Since the running coupling diverges at some infrared 
scale, Aqcd, we have to regularize the running coupling in the low energy region, 
for which there exist several ways (see, e.g., Refs. j2B 1221 123 I2H ESI)- Even if 
we fix the infrared regularization in such a way that we can solve the SD equation 
on the real (space-like) axis, another problem arises when we try to solve the HBS 
equation. Since the argument of the quark mass function in the HBS equations 
for the massive bound states becomes a complex quantity after the Wick rotation 
has been made, we have to solve the SD equation on the complex plane, which 
requires an analytic continuation of the running coupling. Several works such as in 
Refs. j2Hl I2Z| proposed models of running couplings for a general complex variable 
which are consistent with perturbative QCD for large space-like momentum. However, 
they still have branch cuts on the complex plane, and it is a complicated task to obtain 
the solution of the SD equation for a general complex variable. One way to avoide a 
such complication is solving the inhomogeneous BS equation for vertex functions to 
obtain the current correlators in the space-like region which we can fit the mass of 
the relevant bound state to (see, e.g., Refs. [2B1I2H])- Another way might be replacing 
the entire running coupling with an ad hoc analytic function (see, e.g., Ref. 30J ) . 
Anyway, it is impossible to solve the SD equation on the complex plane without 
modeling the running coupling. 

In this paper, we point out that the situation dramatically changes when we 
increase the number of massless quarks. When Nf becomes larger than Nj ~ 8.05, 
the running coupling obtained from the renormalization group equation (RGE) with 
two-loop approximation takes a finite value for all the range of the energy region 
due to the emergence of the non-trivial IR fixed point. Then, we need no infrared 
regularization, and we do not have any ambiguities coming from the regularization 
scheme which do exist in the case of small Nf. Moreover, an explicit solution of the 
two- loop RGE can be written in terms of the Lambert W function jH |M] , and when 
Nf is close to Nf nt the solution of the RGE has no singularity on the complex plane 
except for the time-like axis [31] . Consequently, we can solve the SD equation on the 
complex plane without introducing any models for the running coupling. 

Based on these facts, we solve the HBS equations for the bound states of quark and 
anti-quark in large Nf QCD with the improved ladder approximation in the Landau 
gauge. The mass function for complex arguments needed in the HBS equation is 
obtained by solving the SD equation with the same kernel as that used in the HBS 
equation. We find the solution of the HBS equation in each of the scalar, vector, and 
axial-vector channels, which implies that the scalar, vector, and axial-vector bound 
states are actually formed near the phase transition point. Our results show that 
the masses of the scalar, vector, and axial-vector bound states go to zero as the 
number of quarks Nf approaches to its critical value Nf Tlt where the chiral symmetry 
restoration takes place. This may suggest the existence of a new type of manifestation 
of chiral symmetry restoration in large Nf QCD other than the linear sigma model like 
manifestation and a simple version of the vector manifestation proposed in Ref. |lUj . 

This paper is organized as follows. In section |21 we numerically solve the SD 
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equation with an approximate form of the running coupling, and study the critical 
behavior of the Nambu-Goldstone boson decay constant. In section El we solve the SD 
equation for complex arguments. Section 0] is devoted to summarizing the numerical 
method for solving the HBS equation. Section is the main part of this paper. 
We first solve the HBS equation for the pseudoscalar bound state to show that the 
approximation adopted in the present analysis is consistent with the chiral symmetry. 
We next solve the HBS equation for the scalar, vector, and axial- vector bound states 
to obtain their masses. Finally we give a summary and discussion in section |3 In 
Appendix [X] we solve the HBS equation for the ortho-positronium with a constant 
electron mass to show the validity of the fictitious eigenvalue method. The bispinor 
bases for the bound states are listed in Appendix |B] In Appendix we calculate 
the coupling constants F v , Fa, and G$ of the vector, axial- vector, and scalar bound 
states to the vector current, axial-vector current, and scalar density. We briefly study 
numerical uncertainties in the present analysis in Appendix [DI 



2 Schwinger-Dyson equation in large Nj QCD 

In this section we numerically solve the Schwinger-Dyson (SD) equation for the quark 
propagator with the improved ladder approximation in the Landau gauge, and show 
the critical behaviors of the dynamical mass and the decay constant of the Nambu- 
Goldstone boson. We also show the behavior of the fermion-antifermion pair conden- 
sate (ipip) near the phase transition point. 



2.1 SD equation in the (improved) ladder approximation 

Schwinger-Dyson (SD) equation is a powerful tool to study the dynamical generation 
of the fermion mass directly from QCD (for reviews, see, e.g., Refs. ^3 E2]). The 
SD equation for the full fermion propagator iSp 1 = A{p 2 )f — B{j> 2 ) in the improved 
ladder approximation [221 IS] is given by (see Fig. ^ for a graphical expression) 

(2.1) 

where C% (= N £ N - ) is the second casimir invariant, and g(p,q) is the running 
coupling. The Landau gauge is adopted for the gauge boson propagator. The SD 
equation provides coupled equations for two functions A and B in the full fermion 
propagator Sf- When we use a simple ansatz for the running coupling, g 2 (p,q) = 
g 2 (ma.x(p 2 E , q 2 E )) [221 121], with (p E ,q E ) being the Euclidean momenta, we can carry 
out the angular integration and get A(p 2 ) = 1 in the Landau gauge. Then the SD 
equation becomes a self-consistent equation for the mass function EQo 2 ) = B(p 2 ). 
The resultant asymptotic behavior of the dynamical mass S(p 2 ) is shown to coincide 
with that obtained by the operator product expansion technique [121 HH] • 

However, it was shown in Ref. [TH| that the axial Ward-Takahashi identity is 
violated in the improved ladder approximation unless the gluon momentum is used 
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Figure 1 : A graphical expression of the SD equation in the (improved) ladder approx- 
imation. 



as the argument of the running coupling as g 2 ({j>E — 1e) 2 )- In this choice we cannot 
carry out the angle integration analytically since the running coupling depends on 
the angle factor cos9 = p E ■ <1eI\Pe\\q.e\- Furthermore, we would need to introduce a 
nonlocal gauge fixing JH] to preserve the condition A — 1. 

In Ref. [21], however, it was shown that an angle averaged form g 2 {p 2 E + q%) gives 
a good approximation. Then, in the present analysis we take the argument of the 
running coupling as 

9 2 (pE,q E ) g 2 (pl + q 2 E ). (2.2) 

After applying this angle approximation and carrying out the angular integration, we 
can show (see, e.g., Refs. [321 HHH D3 El) that A always satisfies A{p 2 ) = 1 in the 
Landau gauge. Then the SD equation becomes 

S( x ) = C 2 I dy V S ^ ^ + V) / 2 3) 

where x = p 2 E and y = q\ . Although the choice of arguments in Eq. (|2.2j) 
explicitly breaks the chiral symmetry as mentioned above, it will be shown later that 
the magnitude of the breaking is negligible. 

2.2 Running coupling in large Nf QCD 

In QCD with Nf flavors of massless quarks, the renormalization group equation 

(RGE) for the running coupling a(/x) ( = ) in the two- loop approximation is 
given by 

u,—a(n) = —ba 2 (fx) — ca 3 (/i), (2.4) 

where 

1 1 / iV 2 - 1 \ 

b = - (11N C - 2N f ) , c = — - lON c N f - 3^— N f j . (2.5) 

From the above beta function we can easily see that, when b > and c < 0, i.e., Nf 
takes a value in the range of Nj < Nf < ^N c (N* f ~ 8.05 for N c = 3), the theory 
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Figure 2: Two-loop running coupling (solid line) compared with the approximate form 
in Eq. JUJ (dashed line) for N f = 9. 



is asymptotically free and the beta function has a zero, corresponding to an infrared 
stable fixed point [UI2], at 

b , x 

a, = . 2.6 

c 

Existence of the infrared fixed point implies that the running coupling takes a 
finite value even in the low energy region. Actually, the solution of the two loop RGE 
in Eq. ()2.4|) can be explicitly written [^U EH] in all the energy region as 



at(fj) 



a* 



W(fx ba */eA ba *) + l 



(2.7) 



where W(x) = F 1 (x) with F(x) = xe x is the Lambert W function, and A is a 
renormalization group invariant scale defined by [2| 



A = fi exp 



b a* 



log 



a((j,y 



a(/x) 



b a(/x) 



(2.8) 



We note that, in the present analysis, we fix the value of A to compare the theories 
with a different number of flavors, and that we have no adjustable parameters in 
the running coupling in Eq. (|2.7j) . accordingly (see discussion below). We show an 
example of a(n) for Nf = 9 by the solid line in Fig. El 

The fact that the running coupling is expressed by a certain function as in Eq. ()2.7j) 
implies that, in the case of large Nf QCD, we do not need to introduce any infrared 
regularizations such as the ones adopted in Refs. [2H 122 I2H] for studying real-life 
QCD with small Nf in which the infrared regularization parameter must be chosen in 
such a way that the running coupling in the infrared region becomes larger than the 
critical value a cr = it /A for realizing the dynamical chiral symmetry breaking [21]. The 
running coupling in large Nf QCD takes a certain value in the infrared region for given 
Nf, so that we can definitely determine, within the framework of the SD equation, 
whether or not the dynamical chiral symmetry breaking is realized. Actually, the 
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value of a* decreases monotonically with increasing Nf, and the chiral symmetry 
restores when Nf becomes large enough. In Refs. 0111, it was shown that the phase 
transition occurs at iVj nt ~ 11.9 for N c = 3 (corresponding to a* = a cr = vr/4). 

In order to reduce the task of numerical calculations in solving the HBS equation, 
we modify the shape of the running coupling. Since the dynamics in the infrared 
region governs the chiral symmetry breaking, we adopt the following approximation 
for the running coupling [3J H] : 

~^y± = a*6(A*-(x + y)). (2.9) 

In this approximation the coupling takes the constant value a* (the value at the IR 
fixed point) below the scale A and entirely vanishes in the energy region above this 
scale. The dashed line In Fig. El represents the approximated form of the running 
coupling for Nf = 9. 

2.3 Numerical solution for the SD equation 

In this subsection we briefly explain how we solve the SD equation numerically. 
We first introduce the infrared (IR) cutoff Xsd and ultraviolet (UV) cutoff A$d 

as 

A 2 e A SD /A < X) y < A 2 e A SD /A_ (2.10) 

Then, we discretize the momentum variable x and y into N$d points as 

Xi = A 2 exp [\ SD /A + DsD-t] , (i = 0, 1, 2, • • • , (N SD - 1)) , (2.11) 

where 

n (A S p - A SD )/A 

DSD = N SD -l ■ (2 ' 12) 
Accordingly, the integration over y is replaced with a summation as 

Jdy D SD J2vj ■ (2-13) 

3 

Then, the SD equation in Eq. (|2.3J) with the running coupling in Eq. (|2.9jl is rewritten 

as 

*«) - ^^E^^ss^y^fc- (2 - 14) 

This discretized version of the SD equation is solved by the recursion relation: 

*W)(*0 = -^Dsn^fixi + yj) ; • (2.15) 

4tt 2 y max(xi, y-) Vj + Ef n) {yj) 

Starting from a suitable initial condition (we choose E(o)(xj) = 1), we update the 

mass function by the above recursion relation. Then, we stop the iteration when the 
convergence condition 

„2 



D SDj2j^[Z(n+l){ X i)- E (n){Xi)\* < S 2 A 6 (2.16) 
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Figure 3: A solution of the discretized SD equation in Eq. (|2.14() for Nf = 9. 



is satisfied for sufficiently small e, and regard this E( n ) as a solution of Eq. (|2.14|) . 
In Fig. |2l we show the numerical solution for the mass function E(x). Here, we took 
Nf = 9 (a* ~ 5.2) as an example and used the following parameters: 

A SD /A = +15 , A SD /A = -15 , iV 5D = 1000 , e = 1(T 15 . (2.17) 

Now, let us study the critical behavior of the fermion mass as Nf is varied. Note 
that we can use a* instead of Nf as an input parameter, because once we choose a 
value of Nf, the value of a* is uniquely determined from Eq. (|2.6|) . For example, 
a* = 1 implies Nf = 11.42 and a* = a CT implies Nf = 11.91. We solve Eq. (J2.14)) for 
various values of a* and plot the values of E(m 2 ) in Fig. 0] Here, m represents the 
dynamical mass defied by m = E(m 2 ) . It should be noticed that m is defined in the 
space-like region which does not represent the pole mass of fermion. As we will show 
in section the present fermion propagator does not have any poles and then there 
are no pole masses of fermion. 

We compare this result with the analytic solution j3J E] : 

E(m 2 ) w A exp I n I for a, > a CI . (2.18) 



In the above form there is an ambiguity in the prefactor. Then, we introduce the 
function 

fr(a.) = dAexp|- ? JL= | , (2.19) 



and fit the value of the pre-factor d by minimizing 

^2 | S(x = m 2 ; a,) - h(a,) | 2 , (2.20) 

a* 

in the range of a* G [0.885 : 1]. The resultant best fitted value of d is 

d~4.0. (2.21) 
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Figure 4: Numerical solutions of T,(x — m 2 ) for several values of a* (indicated by O). 
Solid line shows the function in Eq. (|2.19|) with the best fitted value d — 4.0. 



We plot the function in Eq. (|2.19|) with the best fitted value d — 4.0 in Fig. 0] (solid 
line). This clearly shows that the a*-dependence of the resultant S(m 2 ) from our 
numerical calculation is consistent with the analytic result: The dynamical mass 
function vanishes when a* reaches the critical value a CT = tt/4. Noting that de- 
creasing a* corresponds to increasing Nf for fixed N c as we discussed in the pre- 
vious subsection, we see that the chiral symmetry restoration actually occurs at 
Nf = Nf 1 ~ 12 (JV c /3) El- 



2.4 Pseudoscalar meson decay constant in large Nf QCD 

So far we have solved the SD equation and obtained the mass functions for the various 
values of a*. Now we can calculate the pseudoscalar meson decay constant Fp at each 
a* by using the Pagels-Stokar formula PS] : 

N n r . x E 2 (x) - 4- -f\Z 2 (x) 



Fl = — dx 4 dx i - w j / 2 22) 

P 4vr 2 J ( x + S 2 (x) f ' K ' J 

In Fig. we plot the values of F P for several values of a* (indicated by O). To study 
the critical behavior of the pseudoscalar meson decay constant we use the function of 
the form in Eq. (j2.19j) and fit the value of d by minimizing 

J2\F P (a*) - h(a*)\ 2 , (2.23) 

for a* G [0.885 : 1]. The resultant best fitted value of d is 

d ~ 1.5 (= d Fp ) . (2.24) 
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Figure 5: Values of Fp calculated from the Pagels-Stokar formula for several values 
of a* (indicated by O). The dotted line shows the function of the form in Eq. I|2.19[) 
with the best fitted value d — 1.5. 



We plot the fitting function with d = 1.5 in Fig. (dotted line). This shows that the 
results of the numerical calculations for Fp are well fitted by the function of the form 
in Eq. (|2.19p . and that the pseudoscalar meson decay constant has the same critical 
behavior as the mass function has. 

2.5 Fermion-antifermion pair condensate 

In this subsection, we calculate the fermion-antifermion pair condensate (ipip) in large 
Nf QCD, and show that the system in the present analysis has the large anomalous 
dimension 7 m for the operator ipip. We also show that the values of 7 TO are not affected 
so much by the approximation for the running coupling used in the present analysis. 
The condensate (tpip) is calculated from the following equation: 

it ,\ N c f k uv , x E(x) 

where is the mass function obtained from the SD equation and huv represents 
UV cutoff introduced to regularize the UV divergence. In the improved ladder approx- 
imation, the high-energy behavior of the mass function is consistent with that derived 
using the operator product expansion (OPE). The chiral condensate calculated using 
the mass function was shown to obey the renormalization group equation derived with 
the OPE (see, e.g., Refs. [HHEE]). Then, as was adopted in Refs. |2J El H] , we identify 
the condensate, which is calculated with UV cutoff Auv, with that renormalized at 
the scale Ajjv in QCD. 1 



^^When the condensate is calculated using the approximated running coupling defined by Eq. 1)2.9(1 . 
the integration in Eq. H2.25|l is effectively cut off at the scale of A due to the truncation of the running 
coupling for any values of Auv > A. (See Fig. 0: T,(x) = for x > A 2 .) 
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We expect that infrared dynamics in large Nf QCD is similar to that of strong 
coupling QED or walking gauge theories (36J since the running coupling in large Nf 
QCD is well approximated by the constant coupling (see Fig. |2j) |4 a . Then, we also 
expect that the value of the anomalous dimension in large Nf QCD becomes 7 m ~ 1 
since the walking gauge theories have 7 m ~ 1 [SEj- 

When a considering system has the anomalous dimension 7 m , scaling properties of 
Fp and — (ipip) with respect to a* near the critical point are expressed as follows 



F P ~ m, (2.26) 

-(^)~m 3 ~ 7m , (2.27) 

where m represents the dynamical fermion mass. These equations mean that the 
relation between (ipip) an d Fp can be written as 

- (iji/j) = c-Fp" 7m , (2.28) 

where c is a certain positive constant. From this equation, we can express the anoma- 
lous dimension as 

lm = l' m - e, (2.29) 



where 



log(-(V^)) 

7 - = 3 - logFp ' (2 ' 30) 
l0gC (2.31) 



log F P 



Here, we note that 7^ approaches 7 m for a* — > a CT since F P becomes small, i.e., 
(— logFp) 3> 1, near the critical point (see Fig. |3J): 

e — ► for a* ^ a CI . (2.32) 

In Fig. we plot the values of 7^ for several values of a* as an estimation of 
the anomalous dimension. The data indicated by O in Fig. |B] is obtained with the 
approximated running coupling (dashed line in Fig. |2J) in the SD equation. (We call 
this kind of data 7'm PF> M O n the other hand, the data indicated by + is the result 
from the calculation with the two-loop running coupling given in Eq. ([2.7)1 . (We 
call this kind of data 7^°" app ' .) 2 From these results, we conclude that large Nf 

2 The reason why we introduced the approximated running coupling l|2.9|l in the present analysis 
is to reduce the task of numerical calculations in solving the HBS equations. As for the SD equation, 
we can easily solve it numerically with the two-loop running coupling given in Eq. I|2.7|l . 

Since we have to compare 7'^ pp ^ and 7'TO°" app ^ at the same energy scale, we have lowered the 
scale of (^V)^ no ~ app ^ from Aj/y to A by the following two- loop renormalization group equation : 



~a(A uv y 


4ivb 


a(A) 





2t^\ a(A) - a(A uv ) 

4-7T& b 2 ) 47T 
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Figure 6: Values of j' m for several values of a*. O and + represent the values of j' m 
calculated from Eq. i|2.30fl with and without approximation for the running coupling. 
(We call them 7'^ PP ^ an d 7'^ OHlpp ' respectively.) 



QCD with two-loop running coupling as well as with approximated running coupling 
actually possesses 3 

Im^l'm- 1- (2-33) 

Moreover, Fig. |U1 shows that the data of 7 / ^ PP ' ) is i n good agreement with that 
of y^°" app ' ( which implies that the approximation for the running coupling used 
in the present analysis works well. We also expect that the approximation does not 
affect the results so much when we calculate the HBS equations for the bound states. 



3 SD equation on the complex plane 

As we will discuss in section HJ we need the mass function for complex arguments 
when we solve the HBS equation for the massive bound state. In this section, we 
first introduce the SD equation for the complex argument following Ref. [37] (see also 
Ref. [38J), and then solve it in the case of large Nf QCD. 

The SD equation for the complex argument is expressed as jSZj 

£(*)=C 2 -^ 

l07T Z 

where, 

7 W = 6C7 2 , 7 « = C 2 (iC 2 + 9 -l Nc - f Nf ^j . 

3 From the values of c obtained by fitting to the data of ("ipip), we find e = 0.04 ~ 0.06 for the 
approximated running coupling, and e = 0.16 ~ 0.25 for the two- loop running coupling. 



/ dy— + dy 

JC(Q,x) X JC(x,oo) 



C(x,oo) 



g 2 (x + y)Z{y) 
V + Z 2 (y) 



(3.1) 
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where C(a,b) is the integral path on the complex plane. Here, we took the same 
argument of the running coupling as that in Eq. (|2.2jl . and carried out the angle 
integration. Note that the integral path C(a,b) must be taken so as to avoid the 
branch cut appearing in the integral. 

We first study the structure of the running coupling appearing in the SD equation 
(13. lj) to clarify the branch cut. In the improved ladder approximation it is essen- 
tial to use the running coupling determined from the /3-function in the high energy 
(space-like) region for consistency with perturbative QCD. In QCD with small Nf, 
however, the running coupling obtained from the perturbative /3-function diverges at 
some infrared scale, Aqcd- In the ordinary SD equation in the space-like region, the 
infrared singularity is avoided by introducing infrared regularization such as the so- 
called Higashijima-Miransky approximation (2UI22] and its extension as in Ref. |23j . 
However, since the running coupling in Eq. (|3.1|) is a complex function which has the 
complex argument, we need an extension with analyticity satisfied. Several works 
such as in Refs. I21j proposed models of running coupling which are consistent 
with perturbative QCD in the high energy region. But they still have branch cuts 
on the complex plane, and it is a burdensome task to evade all the branch cuts by 
carefully selecting the integral path in Eq. (|3.1|) . One way to avoide such a complica- 
tion might be to give up the consistency with perturbative QCD and use models of 
running coupling with analyticity such as the one used in Ref. |3*U] . 

Here we point out that the situation dramatically changes in the large Nf QCD. In 
the case of large Nf QCD, as we explained in subsection 12.21 the running coupling, as 
well as the two-loop /^-function, is finite for any space-like momentum. This implies 
that we may be able to construct the running coupling by analytic continuation using 
the /3-function. Actually, an explicit solution of the two-loop renormalization group 
equation (RGE) can be written in terms of the Lambert W function |H . When 
Nf is close to A^ nt , the solution of the RGE has no singularity on the complex plane 
except for the time- like axis |3*T| . 

As a result, for general complex x except on the time-like axis (x < 0), we can 
take the integral path C(a, b) in such a way that it just avoids the branch cut coming 
from the angle integration. In Fig. [7| we show the branch cut together with a simple 
choice of the integral path [37|. We stress again that the reason why we can take this 
simple integral path is that the running coupling has no singularity on the complex 
plane except for the time-like axis. 

For solving the SD equation on the complex plane, we here study the explicit form 
of the running coupling. In Fig. |S] we show the real part of the running coupling on 
the complex plane which is obtained by performing the analytic continuation from 
the running coupling on the real axis determined from the two-loop /3-function. This 
figure shows that Rea ~ a*, i.e., Ima ~ 0, in the range of Y — \y\ < A 2 , and that 
Rea <C a* in the range of Y > A 2 . Thus we take the following approximation for the 
running coupling on the complex plane: 

a(y) = a* 6(A 2 - Y) , (3.2) 

which is smoothly connected to the approximation adopted in Eq. (|2.9|) for the running 
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Figure 7: Integral path of the SD equation (|3.1|) . Here, the branch cut appears from 
the four-dimensional angle integration. 




Figure 8: Real part of the two-loop running coupling for Nf — 11 on the complex 
plane obtained by the analytic continuation from the running coupling on the real axis 
(we use the Cauchy-Rieman relation). The complex argument of a is expressed as 
y = Ye l9 , where Y and 9 are real. Note that y is in the space-like region for = and 
in the time-like region for 8 = ir. 
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coupling on the space-like axis. 

Now, let us solve the SD equation 1)3.1)) to obtain the mass function for com- 
plex variable x. Along the integral path shown in Fig. the variables x and y are 
parametrized as 

x = Xe i0 , y = Ye i0 , (3.3) 
where X, Y, and 9 are real. Then the SD equation (j3.1)) is rewritten as 



i0\ 



Co 



16tt 2 



x y 
dY— 
o X 



dY 



x 



i0\ 



f(X + Y) S(Fe 
Y + e~ id T, 2 (Ye i9 ) 



From this we can easily see that the solution is expressed as 

Ux) 



£(X) 



where is real and satisfies the original SD equation on the real axis: 



S(X) = c 2 



16tt 2 



dY—+ / 
o X Jx 



g 2 (X + Y)E(Y) 
Y + E(Y) 2 



(3.4) 



(3.5) 



(3.6) 



Note that the fermion propagator Sf does not have any poles since the kinetic part 



x and the mass part £ (x) have the same phases as x + S z 



x 



e^(X + £ 2 (X)) [see 



Eq. ()3.5)) ] and the mass function in the space-like region satisfies X + Y?(X) > 0. 

We should note that the above solution in Eq. ()3.5j) is a double-valued function 
on the complex plane: The variable x = Xe ld can be parametrized as x = Xe ! ' 9+2?r ' 
for which the solution takes S(x) = e i( - 6, / 2+7r ^S(X) = — e l6l / 2 E(X). This corresponds 
to the fact that the SD equation has two solutions: When S(x) is a solution, — E(ar) 
also satisfies the equation. When we choose the range of 9 as 6 E [—n : ir], the branch 
cut emerges on the time-like axis. This choice is natural because the appearance of 
the branch cut in the time-like region seems consistent with the analytic structure of 
the running coupling. We will see that this branch cut does not matter in calculating 
the bound state masses. 



4 Homogeneous Bethe-Salpeter equation 

In this section we briefly review the homogeneous Bethe-Salpeter (HBS) equation for 
the bound states of quark and antiquark and show how to solve it numerically. 

4.1 Bethe-Salpeter amplitude 

In this paper, we consider the scalar, pseudoscalar, vector, and axial-vector bound 
states of quark and antiquark, and we write these bound states as \S(q)), \P(q)), 
\V(q, e)), and \A(q, e)), respectively. Here, g M represents the momentum of the bound 
states and e M represents the polarization vector satisfying e • q = and e 2 = —1. 
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Now, we introduce the Bethe-Salpeter (BS) amplitudes x f° r the bound states of 
quark and anti-quark as follows : 



(0| T ^ afi (x+) \S a (q)) = 5: 



T^ afl (x + ) jp(x-) \P a (q)) = & 



c -iqX 
.j (Ag)/ _ iqX 



V2 
f (A« 



d A p 
d 4 p 

W) 4 



e- ipr [X(s)( P ;q)U, (4.1) 



-vpr 



lx(p)(p;q)U, (4.2) 



(0| T ip afi (x + ) ffi(z-) \V a (q,e)) = 5{ K -^J-e^ x 

p -iqX 

V2 



d 4 p 
(2^ 

d 4 p 



-ipr 



[X(v)(p] q,e)) a /3, 



(4-3) 



-vpr 



[X(A)(p; q,e)]ap, 



(4-4) 

where x± = X ± r/2, A a is the generator of SU(iVj) normalized as tr[A a Ab] = 25 a b, 
and (a, /?), (/, /'), and (z, j) denote the spinor, flavor, and color indices, respectively. 

We can expand the BS amplitude \ in terms of the bispinor bases P and the 
invariant amplitudes \ l as follows : 



X(s,p) (p; ?)] = J2 [ T \s,p) (p; ?)] , X(5,p) (p; ?) , 



i=i 

8 



X(V,A)(P> ?> 



a/3 



Z [ r V,A) (p; ?> e )] Qj3 xV,A) (p; ?) • 



(4.5) 



(4.6) 



i=i 



The bispinor bases can be determined from spin, parity, and charge conjugation 
properties of the bound states. The explicit forms of T(p), ^\v)i anc ^ ^(A) are 
summarized in Appendix [Bj 

We take the rest frame of the bound state as a frame of reference: 



g" = (M B , 0,0,0), 



(4.7) 



where Mb represents the bound state mass. After the Wick rotation, we parametrize 
p M by the real variables u and x as 

p ■ q = iM B u , p 2 = —u 2 — x 2 . (4.8) 

Consequently, the invariant amplitudes \ l become functions in u and x: 

X\ S ,P,V,A) = x\s,P,VA)( U i X )- ( 4 - 9 ) 

From the charge conjugation properties for the BS amplitude \ and the bispinor 
bases defined in Appendix |E1 the invariant amplitudes x l { u -, x ) are shown to satisfy 
the following relation: 



x\s,p,v,a)(, u ' x ) — x\s,py,A)(~ u i x ) ■ 



(4.10) 
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Figure 9: A graphical representation of the HBS equation in the (improved) ladder 
approximation 

4.2 HBS equation 

The HBS equation is the self-consistent equation for the BS amplitude (see, for a 
review, Ref. [H]), and it is expressed as (see Fig. 

T X = K X ■ (4.11) 

The kinetic part T is given by 

T( P -q) = iS- F \ P + q/2)®iS F \ P -q/2) , (4.12) 

where Sf is the full fermion propagator ( iS F 1 (p) — j) — £), and the BS kernel K in 
the improved ladder approximation is expressed as 

In the above expressions we used the tensor product notation 

(A®B) X = A X B , (4.14) 

and the inner product notation 

K X (p; q)= J K(p, k) X (k; q) . (4.15) 

It should be noticed that the fermion propagators included in T in Eq. ()4.12j) 
have complex-valued arguments after the Wick rotation. The arguments of the mass 
functions appearing in two legs of the BS amplitude are expressed as 



- (p±q/2) 2 = u 2 + x 2 - (^j-J TiuM B . (4.16) 

In general, it is difficult to obtain mass functions for complex arguments. However, 
as we have shown in sectional it is easy to obtain them in the case of large Nf QCD. 
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We now modify Eq. (|4.11|) so that we can solve it numerically. 4 It is convenient 
to introduce the conjugate bispinor bases defined by 



r(p;q,e)= lo r(p*;q,e)h . 



(4.17) 



Multiplying these conjugate bispinor bases from the left, taking the trace of spinor 
indices, and summing over the polarizations, we rewrite Eq. 1)4.11)1 into the following 
form: 

f y'^ 

T ij (u,x)x 3 (u,x) = I — ^-g — K ij ( y u,x;v,y)x J (v,y), (4.18) 



where the summation over the index j is understood, and 



Tij(u,x) 



K i:j (u,x;v,y) 



Yl \ tT [ r *(P5 9, t)T{p; q)T J {p; q, e) 
1 dcosO Y.\tT[r(p-q,e)K(p,k)P(k;q,e) 



(4.19) 
(4.20) 



with the real variables v and y introduced as k ■ q = wMb and k ■ p = —uv — xy cos 9. 
Here, we note that although the mass function S(a;) has the branch cut on the time- 
like axis as mentioned in sectional has no singularity and becomes a continuous 
function for all the range of u and x. As for K^, the branch cut of running coupling 
g does not matter since its argument (p 2 E + k 2 E ) never takes a negative value. 

Using the property of x 1 m Eq. (J4.10)) . we restrict the integration range as v > 0: 



dvK ij (u,x;v,y)x J (v,y) 



v>0 



dv [Kij(u, x; v, y) + Kij{u, x; -v, y)) x J (v, y). 



(4.21) 

Then, all the variables u, x, v, and y can be treated as positive values. 

To discretize the variables u, x, v, and y we introduce new variables U, X, V, and 
Y as 



u = A e v ' k 
v = A e v ' K 



x = A e x ' h 
y = A e y / A 



and set UV and IR cutoffs as 

U,V e [Xu,Au}, X,Y e [X x ,Ax] 



(4.22) 



(4.23) 



We discretize the variables U and V into Nbs,u points evenly, and X and Y into 
Nbs,x points. Then, the original variables are labeled as 

U[i v ] = Aexp [Xu/A + Dulu] , x [Ix] = Aexp [Ax/ A + D X I X ] , 
V[ Iv ] = Aexp [Xu/A + Duly] , y [Iy] = Aexp [Ax/A + D x Iy] , 



4 In the following we explain the method for the vector and axial- vector bound states. The 
extension to the scalar and pseudoscalar bound states are easily done. 
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where Iu, I v = 0, 1, 2, • • • (N B s,u - 1) and I x , h = 0, 1, 2, • • • (N BS ,x - 1)- The mea- 
sures and D x are defined as 



(At, - AcQ/A 
^(7 = — 77 ; — J 



(A x - A x )/A 
Nbs,u — 1 Nbs,x — 1 

As a result, the integration is converted into the summation: 

y 2 dy dv ■■■ =^ n n x ^ — 3 



v>0 



DuD v v v 3 

IvJy 



In order to avoid integrable singularities in the kernel K(u, x; v, y) at (u, x) 
we adopt the following four-splitting prescription 



Kij(u,x,v,y) 



(4.24) 

(4.25) 

0,2/), 



- [ K i:j (u, x, v+, y+) + Kij(u, x, v+, y_) 
+ Kij(u, x, v-, y+) + Kij(u, x, u_, y_) ], 



where 



u ± = exp 



y± = exp 



y ± 



D 



x 



(4.26) 
(4.27) 



4.3 Fictitious eigenvalue method 



Now that all the variables have become discrete and the original integral equation 
(14.11)) turned into a linear algebraic one, we are able to deal it numerically. However, 
it is difficult to find the bound state mass M B and the corresponding BS amplitude 
X directly since the HBS equation depends on Mb nonlinearly. A way which enables 
us to solve the nonlinear eigenvalue problem is the fictitious eigenvalue method |14j . 
In this method we introduce a fictitious eigenvalue A and interpret the HBS equa- 
tion (|4.11|) as a linear eigenvalue equation for a given bound state mass M B : 



T X = A • K X . 



(4.28) 



Consequently, the HBS equation turns into an ordinary eigenvalue problem which we 
can solve by standard algebraic techniques. By adjusting an input mass M B such that 
an eigenvalue A equals unity, we obtain the bound state mass and the corresponding 
BS amplitude as a solution of the original HBS equation (|4.11j) . In Appendix El to 
show the validity of this method, we calculate the mass of the positronium using this 
method. 



5 Numerical Analysis 



In this section we show the results of our numerical analysis. 
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a* 


A 


a* 


A 


0.89 


1.00121 


0.95 


1.00262 


0.90 


1.00205 


0.96 


1.00267 


0.91 


1.00230 


0.97 


1.00273 


0.92 


1.00241 


0.98 


1.00279 


0.93 


1.00249 


0.99 


1.00284 


0.94 


1.00255 


1.00 


1.00290 



Table 1: Fictitious eigenvalues obtained by solving Eq. H4.28J1 for the pseudoscalar 
bound state with zero mass used as an input. 

5.1 Pseudoscalar bound state 

As discussed in subsection I2.1| the approximation to the argument of the running 
coupling in Eq. (j2.2j) breaks the chiral symmetry explicitly [18] . So, before solving 
the HBS equation for the massive bound states, we solve that for the pseudoscalar 
bound state and see how much the chiral symmetry is explicitly broken by this ap- 
proximation. 

The mass of the lowest-lying pseudoscalar bound state should become zero because 
it appears as a Nambu-Goldstone boson when the chiral symmetry is spontaneously 
broken. So, we substitute zero for the bound state mass and check whether the 
fictitious eigenvalue A becomes unity. 

We use the following parameters for the calculations: 

[ X v , A v ] = [ -18.0, ], [ X x , A x ) = [ -8.5, ], (5.1) 

N BS ,u = N BS , X = 30 . (5.2) 

We calculate the fictitious eigenvalues for several values of a* and show them in 
TableHJ We can see that A = 1 is satisfied within 0.3% uncertainty. This implies that 
our calculations actually reproduce the massless Nambu-Goldstone boson within the 
numerical error, and that the effect of explicit chiral symmetry breaking caused by 
the approximation for the running coupling is negligible. 

5.2 Vector, axial-vector, and scalar bound states 

In this subsection we show the results of the numerical calculations for the masses of 
the vector, axial- vector, and scalar bound states. For the UV and IR cutoffs we adjust 
the values of them in such a way that the dominant supports of the integrands of the 
decay constant in Eq. (|C4|) and the normalization condition in Eq. (|C.5|) lie in the 
energy region between the UV and IR cutoffs. As an example, we show the integrands 
of the decay constant and the normalization condition for the vector bound states 
in Appendix [01 From these figures, the dominant supports lie in the lower energy 
region for smaller value of a*. Then, we use the following a*-dependent UV and IR 
cutoffs for the vector and the axial-vector bound states: 

[Ac/, Au] = [ -12.0 + 22.0 x (a*- 1.0) , -1.0 + 35.0 x (a,- 1.0) ], (5.3) 
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[A x , Ax] = [ -5.0 + 22.0 x (a* - 1.0) , -2.0 + 20.0 x (a* - 1.0) ]. (5.4) 



For the scalar bound state, on the other hand, we use the following fixed UV and IR 
cutoffs: 

[ A v , A v ] = [ -18.0, ], [\ x , A x ] = [ -10.0, ]. (5.5) 

Although the integrands of the normalization conditions are shown in Appendix [D] 
only for the vector bound states, we have checked that the dominant supports always 
lie within the energy region between UV and IR cutoffs for all kinds of bound states 
and for all values of a*. As for the numbers of the discretization, we use 

N BS ,u = 20 , N BS , X = 55, (5.6) 

for the vector and the axial-vector bound states, and 

N BS ,u = N BS , X = 30, (5.7) 

for the scalar bound states. In Appendix El we show that these numbers of discretiza- 
tion are large enough for the present analysis. 

We should stress that we actually found a solution for Eq. ()4.28|) reproducing 
A = 1 for all the types of the bound states in the range of a* G [0.885 : 1]. This 
means that there do exist the vector, axial-vector, and scalar bound states near the 
phase transition point in the broken phase. 5 

Now, let us show the critical behavior of the masses of the existing bound states. 
In Fig. HOI we plot all the bound state masses calculated for several values of a* 
together with the pseudoscalar meson masses obtained in the previous subsection. 
This figure shows that the masses of the vector, axial-vector, and scalar bound states 
go to zero simultaneously as the coupling approaches its critical value (or, equivalently, 
N f -> A^ rit ) : 

M s , My, M A -> for a, -> a cr . (5.8) 

Next, to study the critical behavior of Ms, My, and Ma we use the function of the 
form in Eq. (J2.19)) and fit the value of d by minimizing 

Y,\M s ,vXa*) ~ Ha*)\ 2 ■ (5.9) 

The resultant best fitted values of d for the scalar, vector, and axial-vector bound 
states are 

d Ms ^ 6.2 , d Mv ^ 16.5 , d MA ~ 17.2 , (5.10) 

respectively. We also plot (the square of) the ratio of the bound state masses to F P 
for several values of a* in Fig. [TT] The dotted lines plotted together with the data in 
this figure represent the values of the following ratios obtained from Eqs. (|5.10|) and 
fl2~2^j) : 

{Mg/Fpf = 17, (My/Fp) 2 = 121 , (M A /F P ) 2 = 132 . (5.11) 

5 On the other hand, we cannot find any solutions for the HBS equations in the symmetric phase, 
i.e., a* < a CT (or, equivalently, Nf > Nj nt ). This fact seems consistent with the property of the 
conformal phase transition which has no bound states in the symmetric phase [30]. 
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Figure 10: Values of the scalar, pseudoscalar, vector, and axial-vector meson masses 
for several values of a* . 



(M A ,v,s,p/Fp) 



200 



150 



100 



50 - 



O (M A / Fp y 

+ (M V /F P f 

□ (m s /f p ) 2 

X (M P /F p f 



^ O O O O O <> O p L<> o M 0_<S, O^O O^sxr- 



-+■+ + +■+ 4- 



oHHHBBB . Q . QCH - H3€] £ ]HE ] HHBB&QQ& 

1 X X X X X X X X X X X X X X X X X X X X X X X X 



f 0.8 
CXcr 



0.85 



0.9 



0.95 



1 



Figure 11: Values of (M A /F P ) 2 , (M V /F P ) 2 , (M S /F P ) 2 and (M P /F P ) 2 for several 
values of a* (indicated by O, +, □, and x, respectively). Dotted lines represent the 
value s of ( M A / F P ) 2 = 17, (M V /F P ) 2 = 121, and (M S /F P ) 2 = 132 obtained from 
Eqs. (flHUjl and lj23i|l . 
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This figure clearly shows that all the masses of the scalar, vector, and axial-vector 
bound states have the same scaling property as that of F P : 



— ~ constant. (5-12) 
Fp 

We can also say that these masses have the same scaling property as that of E(m 2 ) 
since Fp and E(m 2 ) have the same scaling property. Their ratios are summarized as 
follows: 

£ 2 (m 2 ) : M'i : M'( : M'\ = 1 : 2.4 : 17.0 : 18.5 . (5.13) 

One might think that the vector and axial-vector boundstates decay into a fermion 
and an antifermion since M v > 4X 2 (m 2 ) and M\ > 4S 2 (m 2 ). However, this does 
not happen: As we noticed above Eq. (|2.18|) . S(m 2 ) is not the pole mass but the 
dynamical mass defined in the space-like region. Furthermore, as we have shown 
in section |3J the fermion propagators do not have any poles in the entire complex 
plane including the time-like axis where the pole mass of fermion should be defined. 
Thus the vector and axial-vector boundstates do not decay into a fermion and an 
antifermion. 



6 Conclusion and Discussion 

In this paper we first pointed out that, when we solve the Schwinger-Dyson (SD) 
equation in large Nf QCD, we do not need to introduce any infrared regularizations 
for the running coupling since it takes a finite value for all the range of the energy 
region due to the existence of the IR fixed point. In the case of small Nf, we have 
to regularize the infrared divergence of the running coupling, and different schemes 
of regularizations would give different results. Furthermore, it is difficult to find 
the regularization which makes the analytic structure of the running coupling simple 
enough. On the contrary, the solution of the two-loop RGE in large Nf QCD is 
explicitly written in terms of the Lambert W function, and the running coupling 
does not have any singularities on the complex plane except for the time-like axis 
when Nf is close to N^ nt . This significant feature of the running coupling in large Nf 
QCD enabled us to solve the SD equation on the complex plane. 

Then, we solved the homogeneous Bethe-Salpeter (HBS) equations for the scalar, 
pseudoscalar, vector, and axial-vector bound states of quark and anti-quark in large 
Nf QCD with the improved ladder approximation in the Landau gauge. In the 
quark propagator included in the HBS equation, we used the quark mass function 
obtained from the SD equation with the same approximation, which is needed for the 
consistency with the chiral symmetry. 

We first showed that the HBS equation provides the massless pseudoscalar bound 
state in the broken phase which is identified with the Nambu-Goldstone boson asso- 
ciated with the spontaneous breaking of the chiral symmetry. Next, we showed that 
there actually exist vector, axial- vector, and scalar bound states even near the phase 
transition point in the broken phase, and that their masses decreases as the number 
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of massless quarks Nf increases. At the critical point all the masses go to zero, show- 
ing the same scaling property as that of the pseudoscalar meson decay constant Fp 
consistently with the picture expected from the conformal phase transition |3| I13|. 

Let us discuss the pattern of the chiral symmetry restoration by considering the 
representation of chiral SU(./V/)l x SU(A^/)r of the low-lying mesons extending the 
analyses done in Refs. jHHl HO] ■ 

For Nf = 3 the pseudoscalar meson denoted by 7r and the longitudinal axial- vector 
meson denoted by Ai are an admixture of (8 , 1) © (1 , 8) and (3 , 3*) © (3* , 3), since 
the chiral symmetry is spontaneously broken (HHl EDJ 

|tt) = |(3, 3*) © (3*, 3))sin^ + |(8, 1) © (1, 8))cosV , 
|Ai(A = 0)> = |(3, 3*)©(3*, 3))cos^- |(8, 1)©(1, 8))sin^ , (6.1) 

where A denotes the helicity in the collinear frame, and the experimental value of 
the mixing angle ip is given by approximately ip = 7r/4 On the other hand, the 
longitudinal vector meson denoted by p belongs to pure (8 , 1) © (1 , 8) and the scalar 
meson denoted by S to pure (3 , 3*) © (3* , 3): 

|p(A = 0)> = |(8, 1) e (1 , 8)> , 

\S) = |(3, 3*) ©(3% 3)) . (6.2) 

When the chiral symmetry is restored at the phase transition point, it is natural 
that the chiral representations coincide with the mass eigenstates: The representa- 
tion mixing is dissolved. From Eq. (jfi.lj) one can easily see that there are two ways 
to express the representations in the Wigner phase of the chiral symmetry: The con- 
ventional manifestation a la the linear sigma model (called the GL manifestation in 
Ref. [TT] ) corresponds to the limit ip — ► tt/2 in which it is in the representation of 
pure (N f , Nf) © (Nf , N f ) of SU(N f ) L x SU(JV/) R together with the scalar meson, 
both being the chiral partners: 

rrn / \*),\S) - \(N f , Nf) © (Nf , N f )) , 

[ ' \ |A 1 (A = 0)),|p(A = 0)) - \(Nj- 1,1)0(1,^-1)). {b - 6) 

On the other hand, the vector manifestation (VM) [TD] corresponds to the limit ip — *■ 
in which the A\ goes to a pure (Nf , Nf)® (NJ , Nf), now degenerate with the scalar 
meson in the same representation, but not with p in (Nj — 1 , 1) © (1 , Nj — 1): 

fV M) { ^'1^ = °)) - 1)©(1, AT|-1)) , 

lVMJ \ |A 1 (A = 0)),|5) - \(Nf , Nf) © (Nf , Nf)) . 

Namely, the degenerate massless 7r and (longitudinal) p at the phase transition point 
are the chiral partners in the representation of (Nj — 1 , 1) © (1 , Nj — 1). 

Now, what does our result say on the chiral representation of low-lying mesons? 
As can be seen from Fig. El the resultant values of the masses obtained from the 
HBS equation roughly satisfies the following relation [40J: 

M 2 A + Mp = Ml + Ml, (6.5) 
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for all values of a*. This relation holds independently of the mixing angle ip given in 
Eq. (j6.1J) when the low-lying mesons saturate the chiral algebra shown in Ref. |40j . 
Then, it is reasonable to discuss the chiral representation without worrying about 
the influence of the excited states of the bound states. By using the relation tan xjj = 
My I Ms [ID] and the values of My and M$ obtained from the HBS equation in 
subsection 15.21 the value of the mixing angle ip is roughly determined as 

tan %j) = My I M s ~ 3. (6.6) 

This implies that 7r and the longitudinal A\ are still admixtures of the pure chiral 
representation even at the chiral restoration point: 

|tt) - \(N f ,N* f )®(N* f ,N f ) )sin^ + \{N 2 f - 1,1) ®{l,N 2 f -l) > cos V, 
\Ai) -> \(N f ,N})®(N* f ,N f ) )cos^ - \{N 2 f -l, 1)©(1, N}-1) ) sin V- 

(6.7) 

This may suggest the existence of a new type of manifestation of chiral symmetry 
restoration in large Nf QCD which is neither of the GL manifestation nor the simple 
version of the vector manifestation (VM). 
Several comments are in order 

In Appendix [nj we show the calculations of the coupling constants Fy, Fa, and 
Gs of the vector, axial- vector, and scalar bound states to the vector current, axial- 
vector current, and scalar density. The results shows that they also have the same 
scaling properties as Fp. These results indicate that all the dimension-full quantities 
determined by the infrared dynamics have the same scaling properties, as far as the 
(improved) ladder approximation is concerned. 

Although the masses obtained from the HBS equation satisfy the condition ()6.5|) 
needed for the saturation of the chiral algebra, the couplings F P , F v , and Fa do 
not seem to satisfy the first Weinberg's sum rule [2]: Fp + F\ = F v . We have 
not fully understood what this means for the pattern of chiral symmetry restora- 
tion. Apparently, reducing the numerical uncertainty will help us to reach the final 
understanding. 

In the present analysis we did not include the effect from the four-fermion inter- 
action which is induced in the case of 7 m ~ 1 as was conjectured in strong coupling 
QED [22] and was demonstrated in walking gauge theories (121 EH] • It is not clear 
at this moment whether or not the qualitative results in the present analysis will be 
changed when we include such an effect. 

In the present analysis we stressed that the running coupling in large Nf QCD de- 
termined from a two-loop /3-function is expressed as the Lambert W function which 
enables us to solve the HBS and SD equations with mutual consistency near the 
critical point. Apparently, this Lambert W function can be used as an infrared regu- 
larization to solve the HBS and SD equations with mutual consistency in the case of 
QCD with small Nf. It will be very interesting to study meson masses near the chiral 
phase transition in hot and/or dense QCD by using such an infrared regularization. 
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We think that it is important to clarify which effective field theory (EFT) describes 
the new pattern of the chiral symmetry restoration expressed in Eq. ([6.7)1 . Especially, 
it is very interesting to see how the matching between the EFT and the underlying 
QCD with large Nf can be done to determine bare parameters of the Lagrangian of 
the EFT. 
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Appendices 

A Positronium 

In this appendix, to show the validity of the fictitious eigenvalue method for solving 
the HBS equation explained in subsection 14. 3[ we calculate the mass of the ortho- 
positronium which is the vector bound state of the electron and the positron. The 
same analysis was done in Ref. j3H], and here we follow the analysis. 

In the weak coupling limit the HBS equation for the ortho-positronium can be 
solved analytically, and the energy spectrum takes the following form |44| : 

il4 n) = 2m e - A.l 

where m e is the mass of the electron and the positron, and a = 1/137 is the coupling 
constant of QED. 

We use following parameters in our calculation: 

[ Xu, Au ] = [ -18.5, -2.9 ] , [\ x , A x ] = [ -10.8, 2.2 ] , (A.2) 

N BS ,u = N BS ,x = 28 , m e = 137.0 . (A.3) 

(We used the energy scale which satisfies the relation m e a = 1 following Ref. J5S •) I n 
Fig. H21 we show the resultant values of the fictitious eigenvalue A for several values 
of the input parameter My. Finding the point where the smallest A becomes unity, 
we determine the value of the mass of the ground state My as the solution of the 
original HBS equation: 

M ( y ] = 273.99842. (A.4) 
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Figure 12: Fictitious eigenvalues for the positronium. 



From this the binding energy is calculated as 

E (1) = 0.00158. 



(A.5) 



These values are in good agreement with the values M^ 1 ' = 273.99817 and = 
0.00183 derived from Eq. (jA.ljl . This shows that our numerical method works well 
to obtain the mass of the ground state. 



B Bispinor bases for scalar, pseudoscalar, vector, 
and axial-vector bound states 

In this appendix we show the explicit forms of the bispinor bases for the scalar, 
pseudoscalar, vector, and axial-vector bound states. Here, we use the notation q^ = 
q^/Ms with Mb being the mass of the bound states, and [a, b, c] = a[b, c] + b[c, a] + 
c[a, b]. 

Bispinor base for the scalar bound state (J PC = ++ ) is given by 

rj s) = l, T\ s) =^ Y\ s) = $(p.q), rfe) = ^,fl, (B.l) 
and that for the pseudoscalar bound state (J PC = 0~ + ) is given by 

r (P) = 7s, rj p) =j){p-q) 75, T\ P) = $ 7s, Tj P) = i ty, $} 75 . (B.2) 
Furthermore, for the vector bound state (J PC = 1 ) we use 

rfo = £ T 2 (v) = \\tMp-q), r? v) = ~[^|], r&o = ^.AIl. ( B - 3 ) 

T 5 {v) = (e.p), r 6 iv) =$(6-p), Tl v) = i(p.q)(e.p), I* v) = ^J](e ■ p), 
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and for the axial- vector bound state (J PC = 1" 



r 1 

1 (A) 

r 4 

1 {A) 



t 75, 



r 2 

1 (A) 



^,p1 75 , V\ A) = -[{,§} (p-g) 75 , 



3!^' ! 



75, 



r (.4) = ^ (e • P) (P • ?) 7s 



F (A) = ( e • P) (P • ?) 75 

1 



r 6 



(A) 



pX e • P) 7s, 
M(e-P) (p-ff) 75- 



(B.4) 



C Coupling constants to currents and scalar den- 
sity 

In this subsection we calculate coupling constants Fy, F A , and G$ of the vector, 
axial- vector, and scalar bound states to the vector current, axial- vector current, and 
scalar density. They are defined by 



< | -0(0) 7 "Y#)l%e)> 

o|^(o) Y 7s y V>(o) I Mq^) ) 

(0|#)^(0)|% e )) 



<U F A M A e", 

$ab Gs, 



(C.l) 
(C.2) 
(C.3) 



where A a is the flavor matrix normalized as tr[A a A&] = 25 a b- 

By using the BS amplitudes for the vector bound state, Fy is expressed as 



FyMy 



V2iN c 



If 



oo poo 

o Jo 



(C.4) 



In the above expression, the normalization of the BS amplitudes \ % are determined 
by the following normalization condition 



2M V 5 e€ , = iN c 



d 4 p 



dT(p~ q) 
Xfe Q, e) —7777 — x(P5 9> e ') 



<9M, 



(C.5) 



Here, we notice again that T(p; q) has no singularity although the fermion propagator 
Sf has a branch cut in the time-like region. So the integral in Eq. (|C.5J) is well-defined. 
Once we have obtained My and the corresponding BS amplitudes by solving the HBS 
equation, we can calculate F v from Eq. (|C.4|) . We can also calculate F A and Gs in a 
similar way. In Fig. [T3] (a) we show the values of F v , F A , and Gs together with F P 
for several values of a*. To see the scaling properties, we plot the ratio of F v , F A , 
and Gs to Fp in Fig. H31 (b). This figure shows that Fy, F A , and Gs have the same 
scaling properties as that of Fp. 
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Figure 13: Values of (a) F v , F A , G s , and F P and (b) G S /F P , F V /F P , and F A /F P 
for several values of a*. 




(a) (b) 

Figure 14: Integrands of (a) the decay constant in Eq. (|C.4(1 and (b) the normalization 
condition in Eq. (|C.5|) for a* = 0.885. 



D Uncertainties for numerical calculations 

To solve the HBS equation for the bound states numerically, we introduced the UV 
and IR cutoffs and converted the HBS equation into a linear eigenvalue equation by 
discretizing the integral. As we have discussed in subsection 15.21 we adjust the values 
of the UV and IR cutoffs in such a way that the dominant supports of the integrands 
of the decay constant in Eq. (|C.4J) and the normalization condition in Eq. (|C.5J) lie 
in the energy region between the UV and IR cutoffs. In Figs. E] and we show 
those integrands for a* = 0.885 and 1.0 in the case of the vector bound state. These 
figures show that the dominant supports lie in the lower energy region for smaller 
value of a*, and that the present choices in Eqs. (|5.3jl and (|5.4|) covers the supports. 
For other values of a* used in the present analysis we have checked that the dominant 
supports always lie within the energy region between the UV and IR cutoffs chosen 
as in Eqs. ((Qj) and (|Ojl . 

As for the numbers of the discretization, due to the limitation of the computer 
resources we used N BS ,u = 20 and N B s,x = 55 as shown in Eq. ()5.6|) . Here we study 
the dependences of the mass and the decay constant of the vector bound state on the 
size of discretization. We show the typical values of the mass in Fig. EI] (a) and those 
of the decay constant in Fig. (b) for five choices of the size of the discretization, 
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(a) 



(b) 



Figure 15: Integrands of (a) the decay constant in Eq. (|C.4J| and (b) the normalization 
condition in Eq. (|C.5|1 for a* = 1.0. 
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Figure 16: Typical values of (a) My j k and (b) Fy/A for five choices of the size of 
discretization, (N BS ,u, N BS ,x) = (14,28), (16,32), (18,36), (20,40), and (20,55). 



(N BS ,u,N BS ,x) = (14,28), (16,32), (18,36), (20,40), and (20,55). Figure ESI (a) 
clearly shows that the choice (Nbs,u, Nbs,x) = (20, 55) is large enough to obtain the 
mass of the vector bound state. On the other hand, Fig. EI] (b) shows that there 
are still uncertainties in the value of the decay constant which come from the size of 
the discretization. Apparently, this uncertainty from the size of discretization is the 
dominant part of the numerical uncertainties in the present analysis. 
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